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Abstract: In this paper, we look for an operator that describes the relationship 
between small errors in representation of the bottom topography in a barotropic 
ocean model and the model's solution. The study shows that the model's so- 
lution is very sensitive to topography perturbations in regions where the flow 
is turbulent. On the other hand, the flow exhibits low sensitivity in laminar 
regions. The quantitative measure of sensitivity is influenced essentially by the 
error growing time. At short time scales, the sensitivity exhibits the polyno- 
mial dependence on the error growing time. And in the long time limit, the 
dependence becomes exponential. 
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Sensibilite d'un modele barotrope de I'ocean 
aux perturbations de topographie. 

Resume : Dans ce papier nous nous interessons a Toperateur qui definit le lien 
entre les petites erreurs dans la representation de la topographie du fond d'un 
modele barotrope ct sa solution. L'etude montre que la solution est tres sensible 
aux perturbations de la topographie dans les regions oii le flux est turbulent. 
D'un autre cote, le flux est peu sensible dans les regions ou recoulement est 
laminaire. La mesure quantitative de la sensibilite est surtout influencee par le 
temps de croissance de I'erreur. Sur de courtes echelles de temps, la dependance 
de la sensibilite au temps de croissance de I'erreur est lineaire, mais sur les 
grandes echelles cette dependance devient exponentielle. 

Mots-cles : Sensibilite, Modele barotrope, Topographie de fond. 
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1 Introduction. 

It is well known the solution of a numerical model depends on a number of 
parameters and parametrisations. One of them, bottom topography, attracts 
much attention because it plays an important role in the determining the flow 
field in the ocean (see for example [5],[T]). However, even if the topography 
is well described, it is not evident how to represent it in the model on the 
model's grid because of the limited resolution. It is known for 30 years, that 
requiring the large scale ocean flow to be well represented, one have to smooth 
the topography to get only corresponding large-scale components of relief [6j. 
In this case, the influence of subgrid-scales has to be parameterized. But it is 
not evident how to do that, and how to apply the parametrisation for a given 
model with a given resolution. It is shown in |15| . that different smoothing of 
the topography pattern may significantly change the model's properties. 

One of possible ways to adapt real bathymetry to a particular model is 
to perform a data assimilation procedure with the topography as the control 
parameter. The control parameter in this procedure is supposed to be modified 
to bring the model within an estimated error of the observations. This idea has 
already been applied in [TT] for a linear shallow-water model in a zonal channel. 

In order to proceed to data assimilation with non-stationary solution of a 
nonlinear model, we evaluate first the sensitivity of model to the topography 
variations. In this paper we address especially the most sensitive and the most 
insensitive modes of the solution with respect to little variations of the topog- 
raphy. In particular, we look for modes to which the solution is not sensitive at 
all. The presence of such modes indicates the impossibility to reconstruct the 
bottom topography from observations in unique way. If the solution exhibits 
no sensitivity to some mode, then the topography may be perturbed by adding 
this mode with no flow change. Mathematically speaking, this mode belongs 
to the null space of the sensitivity operator. The dimension of the null space 
determines the number of independent topography variations resulting in the 
same observable flow. The existence of the null space has recently been pointed 
out by [11 for a shallow water-model on the C grid. In this paper we also 
discuss the null space of the operator that describes the model's sensitivity. 

On the other hand, the most sensitive modes will form the sensitive space. 
Any small perturbation of the topography by a function from this space will re- 
sult in a drastic change of the flow. Concerning the data assimilation procedure, 
it is this space that has to be assimilated in the best and in the fastest way. 
The dimension of this space shows the least number of functions participating 
in the cost functional. 

The model used in this paper is a simple barotropic vorticity equation over 
topography. This model is used in two configurations: a square with flat bottom 
and the North Atlantic region with realistic topography. We discuss the steady 
state solution as well as a non-stationary flow. The sensitivity of an ocean 
general circulation model to bottom topography has already been studied by 
adjoint method in lOj, but using a simple well known model allows us to perform 
a complete study in the whole phase space. In particular, a simplest model 
makes it possible to compare the sensitivity to the topography perturbations 
with the sensitivity to perturbations of other model's parameters, like it's initial 
conditions. 
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The paper is organized as follows. The second section describes the model 
and the sensitivity estimates. In the third section we present the numerical 
discretisation of the problem. The fourth and the fifth sections are devoted to 
experiments in the square box and in the North Atlantic respectively. 



2 Sensitivity estimates. 

We consider shallow-water model with the rigid lid assumption 

du du du 1 dp r^^) 

— - /w + 7i— + w— H — — —-au + vAu 

ot ax oy po ox PoHq 

dv dv dv 1 dp r'-y^ 

— + /M + U— H — = — —-av + vAv (1) 

ot ox Oy Po oy poHo 

d{Hu) d{Hv) _ ^ 

dx dy 

where po is the mean density of water and Hq is the characteristic depth of 
the basin. The Coriolis parameter / is supposed to be linear in y coordinate: 
/ = /o + (3y. 

The third equation allows us to introduce the streamfunction ip, such as 

dtp 

Denoting the vorticity hy uo = ^ — ^ we get 

d I difj d 1 dip 
dx H dx dy H dy 

Using this notation we calculate the curl of the first two equations of the sys- 
tem (P). We get 

— + (lu + Jj dtv u + u h V = vIXuj - crcj H -— (4) 

dt dx dy PqHq 

or 

, I ^ + f \ A , ^i^^y) rK\ 

m + ^^^^ -ir^ = - + 1^ 

where J{ip,uj) — ~ Jacobian operator and T(x,y) = 

_dT^ dTy 
dy dx ■ 

The system ^ is considered in the bounded domain VL and is subjected to 
the impermeability and slip boundary conditions: 

V' \dn= 0, 00 Iso- (6) 

Let us suppose the couple ili{x, y, t), uj{x, y, t) is a solution of the system ([2]), 
([5]) with a given topography H — H[x,y). If we perturb the topography by 
some small 5H, we get another solution of the system {ip + dip , w + 5uj} . 



INRIA 



Sensitivity of a Barotropic Model to the Bottom Topography. 



5 



Our purpose is to define the relationship between 6H and 6uj supposing both 
of them to be sufficiently small: 

\\6H\\ « ||iJ|| and ||<5c.I| « 

We start from the stationary equation So far, the perturbed couple 

{-0 + Sip, uj + (5a;} is the solution of the system with perturbed topography, it 
must also satisfy the equation jS]) 

d 1 dilj + 5i}) 



ui + 5uj — 

dxH + 5H 

Using the Taylor development 

1 1 
H 



d 



1 



9-0 + 5^1) 



dx 



1 - 



5H_ 
IT 



dy H + 5H dy 



H + SH 

and neglecting high order terms, we get from ([7]) 

d ( 1 m\ + ^ d (I 



(7) 



(8) 



5ljj 



dx\H J dx 

The difference between this equation and the equation 



SH\dip + 
dy\H IPJ dy 



5i! 



IS 





--( 


^6H^ 


dSij} 


/ dx 


dy ^ 







So far, both 6ip and SH are supposed to be small, we neglect their product and 
write briefly 

1 5H 
y ySi;^SLJ + V—^Vi} (9) 

This equation allows us to find the perturbation of the streamfunction from 
perturbations of vorticity and topography. 

To get the vorticity perturbation, we consider the evolution equation ([Sjl. 
As well as above, we write the equation for the perturbed topography using the 
development ^ and neglecting high order terms: 



d{uj + Suj) 
di 



> + f Suj 

H 

H H 



H H 



~ vAlu + vASlu — auj — aSuj + curl 



PoHq 



(10) 



The difference between the perturbed equation and the non-perturbed one ([5]) 
writes 

dSuj Suj uj + fo + PySH uj + fo + , r,fx2\ 

-^+^(Va + ) = I'ASuj-crSuj + OiS ) 

in this equation we have also dropped out the term J'{5ip, ^ _ ^ + 'ff 
that contains the product of functions supposed to be small. Using the expres- 
sion for Sijj from ([H]), we get 



dSuj 
"dt 



Suj uj + fa + f3ySH^_^^uj + fo + Py 



H H 



H 



H 



vlv 



Suj) 



H 



S H 

V-^ VV' ) ) = vASuj - aSuj + 0{S'^) 
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This equation can be written in a short niatricial form 

^^A{i^,uj)6cu + B{i^,uj)^-^ (11) 
where operators A and B are defined as 

The system ([TT|) starts from the zero initial state 5[j(a;,y,0) = because the 
purpose is confined on the study of the sensitivity of the solution to the to- 
pography, rather than to the initial state. We require the perturbed solution 
{-0 + Sip^Lo + 5lu} to have the same boundary condition as {ip^uj} because we 
do not want to study the model's sensitivity to boundary conditions. Hence, 
perturbations 5ip,5uj in equations (O, (fTTj) must satisfy 

5^1; \an^ 0, 5uj \an= (14) 

Analysing the form of the equation (fTTj) . we can see the right-hand-side is 
composed by two terms A and B ( (fT2|) . (fT3|) ). The first one, A, is responsible 
for the evolution of a small perturbation by the model's dynamics, while the 
second one, B, determines the way how the uncertainty is introduced into the 
model. The first term is similar for any sensitivity analysis, while the second 
one is specific to the particular variable under study. This term is absent when 
the sensitivity to initial point is studied because the uncertainty is introduced 
only once, at the beginning of the model integration. But, when the uncertainty 
is presented in the bottom topography or in some other internal parameter of 
the model, the perturbation is introduced at each time step and requires an 
additional operator. 

If the reference trajectory is a stationary point, then operators A and B do 
not depend on time. In this case the system ([TT|) can be resolved explicitly. The 
solution is 

5u.[T) -if^_I)£^^-§ = G{T)^-§ (15) 

G(T) 

If the reference trajectory {-ip^uj} is not stationary, discrete time stepping 
is used to integrate (jlip using some numerical scheme. For the simplest Euler 
scheme, iterations procedure on the part of the reference trajectory w(to)...w(io+ 
T) can be calculated in the following way: 

— — = A{^{U),uj{U,))5iu'' + B{^{U),u:{U)Y— 

T H 

Let us suppose that there exists a matrix G„ such as at the nth time step 
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Then, at the n + 1th time step we get 

m _ 

H 



I + TA{ij{tn),i0{tn)) ] Gn + TB{i,{tn),w{tn)) 

5H 



Thus, the matrix G„+i is calculated as 

Gn+l = (l + TA{'li^{tn),U:{tn))]Gn + TB{'ll,{tn),u{tn)) (16) 



for n = 0, 1, ..A'' = T/r with Go = 0. It has to be noted here, the operator Gn 
for a non stationary trajectory depends not only on the time interval T, but also 
on the trajectory part passed by the reference solution 'ip{t),uj{t),t = to . . . to+T. 
So far, the solution of the reference system is imiquc, wo can say the operator Gn 
depends on the initial point ip{to),ijj{to)- To express this dependence explicitly 
we shall write the operator Gn as G{to,T). 

In this paper we arc looking for the most dangerous pcrtur})ation of the 
topography. That means, the topography's perturbation of a given small norm 
which maximizes the norm of the solution's perturbation \\6ui\\ at the time T. 
We can chose any norm in this consideration 

WSu)]]"^ =^)CSuj{x,y,t),5uj{x,y,t)» / ICLj{x,y,t)ui{x,y,t)dxdy (17) 

Jq 

with an auto-adjoint positive definite operator IC. Thus, for example, the energy 
of the solution is obtained with the inverse Laplace operator /C = A~^, the 
enstrophy corresponds to identity operator JC = I. 
In other words we are looking for the 

\\SLu{T)f «ICSuj{T).Slu{T)» 

= '""^ «IC6H,6H» = 

«G*{to,T)ICG{to,T)SH, 5H» 



= max 



max 



«ICSH,SH» 
«KlC-^G*{to, T)ICG{to, T)SH, SH» 



«ICSH, SH» 
= maxX'^{IC-^G*{ta,T)ICG{to,T)) (18) 

where X^{G*{to,T)G{to,T)) arc defined as eigenvalues of the problem 

IC-^G*{to,T)ICG{to,T)ipi = A^^i 

Below the enstrophy of the solution is considered as its norm 

\\6Luf= I Lo^{x,y,t)dxdy (19) 
Jn 

and the corresponding eigenvalues problem 

G*{to,T)G{to,T)<pi = X^i^i (20) 
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The eigenvalues show the growth rate of different modes of the perturba- 
tion. Taking into account the expression for G{to,T) in the stationary case 
()15p we see that in the infinite time hmit T — > oo, the operator tends to 
B*(A*)-i(e'^^ - /)*(e^^ - I)A-'^B. The only dependence on T is associated 

with the exponential. If we consider the expression lirriT >oo^^t^ we see the 

value of this expression tends to the liniT roc^^j^, where are eigenvalues of 

the pure exponential operator e^^: 

That means the growth rate of the perturbation is determined by the maximal 
Lyapunov exponent of the model, i.e. it is the same rate as in the study of 
sensitivity to initial state. This reasonable conclusion shows the system develops 
its own instability. No matter what was the source of perturbation, it will behave 
according to the system's internal instability modes on long time scales. 

The same behavior can be expected in the non stationary case. Despite there 
is no explicit exponent in the expression p6p . it is well known the nonlinear 
model on a strange attractor reveals exponential growth rate of perturbations 
on infinite time scales. So, the intrinsic instability of the non stationary model 
will also dominate on long times. 

However, on short time scales the instability will show different behavior. 
This is natural because the matrices and B in (fTS]) are comparable with the 
exponent on short time scales, i.e. the source of the perturbation is important 
on these scales. 



3 Time dependent solution 

The model has been discretised in space using finite elements method. Details 
of discretisation and construction of matrices A and B ( ([T^ . (fT5|) ) are shown 
in the Appendix. 

The grids used in this paper are presented in figHJ The triangulation of the 
square is composed of 206 triangles. The integration points set, being a union of 
vertices and mi-edges of triangles, counts 445 nodes. The resolution of the grid 
varies between 1/80 of the side length (about 50 km for the square of 4000 km) 
near the western boundary and 1/10 of the side length (about 400 km) near the 
eastern one. 

The triangulation of the North Atlantic is also performed with a grid refine- 
ment near the American coast. This triangulation is composed of 195 triangles 
and 436 points. The resolution of this grid is about the same, i.e. about 40 km 
near the American coast and about 400 km near the European one. 
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Figure lA. Triangulation of the Figure IB. Triangulation of the North 

square. Atlantic region. 

In the experiment with the square box we take the characteristic length of 
the basin L = 4000 km. The bottom is supposed to be either flat or with 
regular sinusoidal topography. We take a steady zonal wind as forcing with now 
classical two gyre antisymmetric pattern. This is seen as a schematic pattern 
for the mean curl of the wind stress over the North Atlantic ocean in middle 
latitudes. Its magnitude is equal to 



y) = ]^ sm — (21) 



where rn — 1.1 "^^"^f is the characteristic wind tension on the surface, 
cm 

The coefficient of Eckman dissipation we chose ascr = 5xl0" 
corresponds to the damping time-scale To- = 2 x lO^s ~ 200 days. The lateral 
friction coefficient v has been chosen in order to avoid numerical instability 
which occurs due to the concentration of variability of the model at grid scales. 



2 

This value has been taken tohe v — 500-^, that corresponds to the damping 
time scale = 3 days for a wave of 100 km length. 

Despite this simplified geometry accompanied by now classical test forcing 
([2T|). this "academic" case has been intensively investigated over the last 15 
years in order to study the role of mesoscale eddies in the ocean circulation. 
This configuration helps us to see the sensitivity in a very well described case. 

The second experiment is carried out in a more realistic configuration. The 
domain was chosen to approximate the North Atlantic region. We assume that 
the domain is comprised in the rectangular between T&'^W . . . 2PW in longitude 
and 15" . . . 65''A in latitude. The boundary of the basin corresponds to the 1 
km depth isobath of the ocean. 

To obtain the forcing in this experiment we have used the data set "Monthly 
Global Ocean Wind Stress Components" prepared and maintained by the Data 
Support Section, Scientific Computing Division, National Center for Atmo- 
spheric Research. These data have been prepared by the routine described in 
[1] . From this data set we choose the mean January wind stress components 
and Ty over the North Atlantic based on 1870-1976 surface observations. These 



data are presented on the 2" x 2° grid. 
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The forcing in this experiment is calculated from these data as 



ip = 20° + yx 50"/i,A = 



dTx 1 dry 

dip cos ip dX ' 



(22) 



-40" 



X X 50° /L 

cos ip 



where L = 5500km is a characteristic length of the basin. The spatial configu- 
ration of the forcing is presented in figHjB. 

The bottom topography has been interpolated from the ET0P05 5-minute 
gridded elevation data T^. It is shown in figHK. 

We chose the coefficient of Eckman dissipation to be the same as in the 
previous experiment, cr = 5 x lO^^s^^. The lateral friction coefficient A has 
been chosen in order to avoid numerical instability which occurs due to the 
concentration of variability of the model at grid scales. This value has been 

2 

taken to be A = 300^^, that corresponds to the damping time scale Ta — Q 
days for a wave of 100 km length. 





Figure 2A. Bottom topography. Con- figure 2B. Forcing of the modeL 

tours from 1000 to 6000m, contour interval Contours from -1.5 to 2.4 a interval 
500 m. 0.3. 

We intend to perform several experiments to see the sensitivity of the model's 
solution to the bottom topography. First, we study the dependence of the largest 
singular values of the matrix G{to, T) (fT5|) on the initial point of the reference 
trajectory tg- 

To obtain the initial point in both experiments, the model has been inte- 
grated during 20 years from the zero state. We suppose that after this period 
the spin-up phase is over and the solution of the model reaches its attractor. 
After the spin-up, the model is integrated during 204.8 days with the time step 
0.1 day. This part of trajectory composed with 2048 samples is used as y, t) 
in (fTTjl to construct the matrix G{to,T) following the procedure pE)) . Average 
streamfunction's plots for the square box and for the North Atlantic region are 
shown in figlSK, figdlA. respectively. To see the length of the part of trajec- 
tory under consideration, we plot the kinetic energy versus the enstrophy of the 
solution in figOP, figHp. 

Ek{t) = / 'ip{x,y,t)ujix,y,t)dxdy ^ H{x,y){u'^{x,y,t) + v'^{x,y,t))d^ 
Jn Jn 

rj{t) = [ ij^ix,y,t)dxdy (24) 
Jn 
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Figure 4B. Enstrophy-energy plot for 

Figure 4A. Mean streamfunction in 204.8 days trajectory in the North At- 
the North Atlantic. lantic. 

In this sequence of 2048 samples we choose the set of initial points to spaced 
by T. Each interval between two adjacent points is used to construct the matrix 
G{to,T). Experiments with 2 values of T were carried out: T = 0.8 and 12.8 
days, i.e. the whole sequence was divided respectively to 256 and 16 subintervals 
of 8 and 128 points each. The value of T we shall call the error growing time, 
because it is during this time the perturbation Suj is allowed to grow. Evolution 
of the 5 largest singular values for each subinterval are shown in figlSl figEl as 
functions of the number of the subinterval in the sequence, or in other words, 
of the reference model's time. 
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Figure 5A. Evolution of the 5 largest 
singular values of the matrix G(to, 0.8) for 
< to < 204.8 days in the square box. 

Singular values x 10'^ 



Figure 5B. Evolution of the 5 largest 
singular values of the matrix G(to,l2.8) 
for < to < 204.8 days in the square box. 




Figure 6A. Evolution of the 5 largest Figure 6B. Evolution of the 5 largest 

eigenvalues of the matrix G(to, 0.8) for < eigenvalues of the matrix G(to,12.8) for 
to < 204.8 days in the North Atlantic. < to < 204.8 days in the North Atlantic. 

We see the sensitivity of the model's solution to the topography is almost 
uniform in the square box. The maximum of the first singular value is approxi- 
mately only twice its minimum. This is not the case in the experiment with the 
North Atlantic, where variations have much larger amplitude. We can see in 
figinithere exists a situation with very low sensitivity to the bottom topography. 
Approximately at 125th day all largest singular values in experiments for both 
T — 0.8 and T = 12.8 days have a clear minimum. Comparing the most and 
the least sensitive situations during the same run, we see the maximum of the 
first singular value is 7-8 times higher than its minimum. 

Together with largest singular values of the operator G{to,T) we ana- 
lyze corresponding singular vectors. These vectors represent the most sensitive 
modes of the solution. We calculate these vectors for both error growing times 
T = 0.8 and 12.8 days. 

An example of the most sensitive singular mode for T = 0.8 days and for 
T = 12.8 days in the square box is shown in figEl Only a part of the total 
region, corresponding to the jet-stream near the eastern boundary is presented 
in these figures. One can see that for short error growing time (i.e. small values 
of T) the mode is concentrated at just several points on the grid. However, for 
long error growing time of twelve days, the eigenmode occupies more important 
region. But, in both cases, the eigenmode is concentrated near the jet-stream 
part of the domain. 
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Figure 7A. The most sensitive singular Figure 7B. The most sensitive singular 

mode for the error growing time T = 0.8 mode for the error growing time T = 12.8 
days. Singular value 7.8 X 10~* days. Singular value 9.2 X 10~^ 

On the other hand, there exists also insensitive modes. Their corresponding 
singular values are equal to and any perturbation of the model's topography 
by one of these modes has no impact on the flow. These modes form the kernel 
of the operator. 

One of these modes can be easily seen from a simple analysis of the model 
If we add to the topography H some perturbation which is proportional 
to H itself SH = aH, the model remains the same. Only the third equation 
is multiplied by 1 + a in this case and that does not disturb the equality to 0. 
Hence, the model exhibits no sensitivity to the perturbation SH = aH and this 
mode belongs to the kernel of the operator G{to, T) (fT6|) . 

Another set of insensitive modes results from the boundary conditions im- 
posed on the vorticity uj in the equation ([5]). As topography H and it's per- 
turbation 5H are defined in closed the domain including its boundary. But 
the boundary conditions on oj require that 5lj is equal to zero on the bound- 
ary because both original equation ([5|) and perturbed one (fTO|) must follow the 
same boundary conditions. The discretised operator G(to,T) is represented, 
hence, by a rectangular matrix with Nq strings and N columns, where A'o is 
the number of internal points of the domain f2 and N is the total number of 
discretisation points including boundary. The operator G* {to,T)G{to,T) (^0]) 
is a square N x N operator possessing as many zero eigenvalues as — Nq. 
These kernel's modes are concentrated on the boundary of the domain. 

Considering internal part of the domain, we can also see the difference in 
sensitivity. Least sensitive modes, corresponding to smallest non null singular 
values represent perturbations of topography to which the model is very little 
sensitive. 

This fact is illustrated in figE] where all non null eigenvalues of the operator 
G* {T)G{T) are shown. The beginning of the spectrum is shown on the zoom. 
We can see singular values decrease rapidly for T — 12.8 days. The 8th singular 
value is already 10 times lower than the first one in the square box. In the North 
Atlantic there exists an outstanding very sensitive first mode which singular 
value is 10 times the second one. So, we need just several modes to approximate 
the whole error behavior. For shorter error growing times (like T — 0.8 days) 
this initial decrease is less sharp and the sensitivity is more uniform. 
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In the middle of the spectrum we find a slowly decreasing sequence of singular 
values. This part represents singular modes with low sensitivity. And at the end 
of the spectrum we see the sharp decrease and several modes with outstandingly 
low sensitivity. 



Singular Values ^ .Singular Values 




100 200 300 400 SO 100 150 200 250 300 350 



Figure 8A. Singular values of the ma- Figure 8B. Singular values of the ma- 

trix G{T) for T = 0.1 and T = 12.8 days trix G{T) for T = 0.1 and T = 12.8 days 
in the square box. in the North Atlantic. 

These modes differ from modes in the kernel of the operator. They have no 
intuitive explanation and they are not concentrated on the boundary. In fact, 
they are not localized in space. They occupy almost the whole domain, espe- 
cially regions where the flow is smooth and laminar. However, these patterns 
can not be considered as a grid noise. Their patterns for both T = 0.8 and 
12.8 days are shown in figlH Contrary to figO whole domain is plotted in these 
figures. 




Figure 9A. The most insensitive sin- Figure 9B. The most insensitive sin- 

gular mode for the error growing time gular mode for the error growing time 
T = 0.8 days. Singular value 3.2 X lO-^*^ T = 12.8 days. Singular value 4.9 X lO"!"* 

Similar effects can be seen in the experiment with the North Atlantic. The 
average most unstable eigenmodes for T = 0.8 days and for T = 12.8 days are 
shown in fig llOl One can see, that for short error growing time T = 0.8 the 
mode is also concentrated in a very small region represented by several points 
on the grid. However, for long error growing time of twelve days, the singular 
mode occupies a more important region also. 
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Figure lOA. The most sensitive singu- 
lar mode in the NA for the error grow- 
ing time T = 0.8 days. Singular value 
1.6 X 10"* 
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Figure lOB. The most sensitive singu- 
lar mode in the NA for the error grow- 
ing time T = 12.8 days. Singular value 
3.3 X lO-'""' 























- -/■ 


































70W 60*' bU* JdlV '.lull iiUII 1 





Figure llA. The most insensitive sin- Figure IIB. The most insensitive sin- 
gular mode in the NA for the error grow- gular mode in the NA for the error grow- 
ing time T = 0.8 days. Singular value ing time T = 12.8 days. Singular value 
1.6 X IQ-is 6.5 X 10-15 



3.1 Stationary point 

When the solution is stationary or supposed to be stationary, we can apply the 
formulae (|15p and avoid discrete time integration. In this case, the computa- 
tional procedure becomes much easier. Consideration of the stationary point 
helps us to see the dependence of singular values on particular parameters, 
rather than on variations of the basic trajectory. 

To get a stationary solution of the model, we increase the dissipation coeffi- 

2 

cient ly up to 3000^^. The spin- up time required to reach the stationary point 
is approximately equal to 600 days. We use 800 days spin-up to ensure the 
stationary behavior of the model after the spin-up. The streamfunction of this 
stationary point is shown in fig ll2l One can see perfect antisymmetric pattern 
produced by antisymmetric wind stress. 
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Figure 12. Streamfunction pattern of the stationary point 



We consider first the dependence of the singular values on the error growing 
time T. As it has been noted above, on large time scales, the dependence 
must be close to the exponential, but on short time scales, when the source 
of perturbation is important, the dependence may differ from an exponential. 
To determine the form of dependence and to identify the tiem scale separating 
"long" and "short" times, we calculate the eigenvalues of G* {T)G{T) for the 
operator G 



G{T) 



TA 



I)A-'B 



WibsKceajaperators A and B are define d ffl'" 1^3 "'' GSl) 

Sinfi.yal.#..l . le+02 




le+OO le+01 le+02 



le+00 le+01 le+02 



Figure 13A. 5 singular values of the 
matrix G(T) in the square box versus error 
growing time T for 10~^ < T < 100 days. 



Figure 13B. 5 singular values of the 
matrix G(T) in the North Atlantic versus 
error growing time T for 10^'^ < T < 100 
days. 

In fig[T3]we can see the behavior of 5 different singular values of G as the 
error growing time T increases. The first (largest) value is plotted together 
with the 10th, the 20th, the 50th, the 100th, the 150th and the 200th. One can 
see that the critical error growing time for all values is close to three days in 
both experiments. For lower error growing times the increasing rate is close to 
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polynomial because the rate is linear in logarithmic coordinates 



lnA = alnT + 6i.e. A - 5T' 



la 



Value of a is equal to 1 in both experiments. Hence, the value of A depends 
linearly on T when T is lower than the critical value. The growth of different 
parts of the spectrum of G is uniform. The lines of the growth of different 
singular values are parallel to each other. 

But, when error growing time exceeds three days, the growth rate of singular 
values accelerates and approaches to the exponential. The growth rate of these 
exponentials is very different in different parts of spectrum. In fact, the first 
singular value increases very rapidly, while the growth of all others is relatively 
slow. This leads to the fact that the ratio of the largest to the smallest singular 
value increases as well as the small number of the most sensitive singular modes 
becomes dominant in the whole sensitivity of the model as we have already seen 
in figia 

Thus, only a few singular values and vectors are important in the sensitivity 
analysis on long time scales. 

We have seen that the sensitivity is linear in time on time scales lower than 
3 days. But, on longer time scales, error growth becomes exponential. That 
means, during 2-3 days of the model's time, the perturbation is introduced into 
the model and, after that, it follows the model's dynamics. During the phase 
of introduction, linear transition of perturbation from topography to model's 
variable dominates, resulting in linear dependence on time. But after 2-3 days, 
it is the model's dynamics that governs the error evolution. Being non-linear 
and intrinsicly unstable, the dynamics ensures exponential error growth of a 
perturbation. On these time scales, topography perturbation evolves like any 
other perturbation from any other source. 

One of interesting questions is to find how sensitive would be the model with 
the same parameters but different topography. In particular, it is interesting to 
see how the sensitivity changes when the topography is not flat but represented 
by mountains of different height or different space scales. 

We perform several experiments with the same model parameters as above, 
but with the bottom topography taken as 



The basic state -0,0; used to construct the operator G is obtained as the sta- 
tionary point reached by the model after the 800 days spin-up. 

In the first experiment, the mountains height a is allowed to vary. We use 
30 equally spaced values from —300 to 300 m: Q!„ = —300 + 20 x n. Values of 
kx and ky were taken both to be 4. The topography is shown in fig|14B. The 
perfect sinusoidal form of the topography is somewhat disturbed near the top 
and the bottom sides of the square due to very low resolution at those places 
(see figHK). 

In this experiment we get the dependence of five singular values of G on the 
bottom mountains height a ffig |14l A_'). Error growing time T in this experiment 
has been chosen as T = 1 day. 



H{x,y) — 500 m 4- asin{kxTTx) sin^kyTry) 



(25) 



RR n° 6717 



18 



E. Kazantsev 




Figure 14A. 5 singular values of the 
matrix G(l day) as functions of mountains 



Figure 14B. Topography in the exper- 
iment with variable a 



height. 

One can see the most sensitive circulation is observed when the bottom is 
flat. If there are some mountains on the bottom, then the sensitivity is lower. 
However, the curves in figtHK are not symmetric with respect to zero. When a 
is negative the largest singular value is smaller than with the same but positive 
a. 

This difference in behavior can be explained by the streamfunction pattern 
of the stationary solution obtained for particular topography figUni When a is 
positive, topography configuration near the left boundary anti-correlates with 
the streamfunction. That means the vorticity pattern is "in phase" with the 
topography. The positive gyre of the streamfunction is situated over the gap on 
the bottom and amplificated, the negative gyre is over the bottom mountain, 
and is amplificated also. 

When a is negative, streamfunction pattern "correlates" with the topogra- 
phy. The positive anomaly in streamfunction is situated over the hill on the 
bottom, and the negative anomaly over the gap. This results in lower strength 
of gyres. Only a little anomaly remains at the usual position in the center of the 
box, the major part is displaced to the North and to the South, ans is situated 
at more favorable topographic position. 
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Figure 15A. Basic state streamfunc- Figure 15B. Basic state streamfunc- 

tion for a = —300m tion for a = -|-300m 

To explain the difference in the sensitivity of the model we may suppose 
that the sensitivity to the topography is related to the general stability of the 
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model's solution. So far, we consider a stationary point of the model, we can 
evaluate the general stability of the point by the smallest dissipation coefficient 
J/ in ([1]) necessary for the stationary point to exist. It is evident when the 
dissipation is strong, the model has a stationary point as a global attractor. Any 
solution starting from any point will tend to the attractor becoming stationary. 
When the parameter v becomes lower, at some value the solution becomes non- 
stationary. It is this value that we compare with the largest singular value of 
the matrix G in figHBl One can see a clear relationship between the stability 
and the sensitivity of the solution. 
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Figure 16. The smallest u of the stationary point vs largest 
singular value of the matrix G 

In the second experiment we look at the sensitivity of the model with the 
sinusoidal topography of different wavenumbers. We use the same formula (|25p 
for H{x, y) with constant a = 100m but with varying wavenumbers and 
ky from (flat bottom) up to 30 (short scale mountains, at the limit of grid 
resolution). Topography pattern for k^ = ky = 10 and k^ = ky = 30 are shown 
in fig ll7l One can see the topography with the wavenumber 10 is well resolved 
in the whole center of the domain, while the wavenumber 30 is only possible to 
resolve in the small region near the jet-stream. 
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Figure 17A. Topography pattern with Figure 17B. Topography pattern with 

/c^j — — 10 — ^y — 30 

Streamfunction patterns of stationary points for these two topographies are 
presented in figUHl We can see both patterns are no longer antisymmetric as 
seen in fig ll2l However, streamfunction's deformation in this experiment is not 
as drastic as in the experiment with the amphtude variations fig llSI The jet- 
stream is concentrated at the usual place and its intensity remains almost the 
same. 




Figure 18A. Basic state streamfunc- Figure 18B. Basic state streamfunc- 

tion for kx = ky = 10 tion for kx = ky = 30 

The dependence of singular values on the topography's wavenumber is pre- 
sented in fig ll9l One can see that variations of singular values are smaller than 
in the experiment with the amplitude change. In addition, there is no tendency 
when wavenumbers increase. The behavior of singular values exhibits irregu- 
lar variations about values of sensitivity of the model with flat bottom. The 
amplitude of variation is relatively small, remaining in frames that the maxi- 
mal value is less than two times the minimal one. In fact, these variations are 
similar to natural variations of the sensitivity due to temporal evolution of the 
circulation with flat bottom. We can not distinguish a particular tendency due 
to the rugosity of the topography. 
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Figure 19. Evolution of the 5 singular values of the matrix 
G(T) for the model with the topography with different wavenum- 
bers 



4 Conclusion 

We have considered the sensitivity of the vorticity field to the topography per- 
turbations in frames of the barotropic ocean model. Both stationary and non- 
stationary solutions of nonlinear model have been viewed. 

We distinguish the analysis of the quantitative measure of the sensitivity, 
expressed in singular values of the operator, and the qualitative pattern of the 
singular function, that precise the geographical region of particularly sensitive 
or insensitive solution. 

The quantitative measure is influenced essentially by error growing time. 
Longer is the time period during which we allow the error to grow, greater 
the sensitivity is. This conclusions is consistent with numerous studies of the 
model's sensitivity to other parameters. Thus, the predictability studies that 
analyze the sensitivity to initial conditions, reveal the exponential (or close to) 
growth rate. Regarding to this, one can cite [2], [8], [12], [13], [7] and many 
others. In this paper we show, that at short time scales, the sensitivity to 
topography differs from the sensitivity to initial conditions. But, in the long time 
limit, the sensitivity of the solution is the same to any source of perturbation. 
The intrinsic model's instability dominates at these time scales and the source 
of the perturbation is no longer important. 

We have analyzed patterns of the most sensitive modes of the solution. These 
patterns point out the regions where the solution more sensitive to the topog- 
raphy perturbations. The sensitivity is important in regions where the flow is 
turbulent. On the other hand, in regions where the flow is laminar, the solution 
exhibits lower sensitivity to the topography. The conclusion is in agreement 
with the result of [11], where it was shown that the sensitivity is largest where 
current speeds are high. 

Barotropic model in the North Atlantic develops turbulent flow near Ameri- 
can coast and laminar flow near the European one. The sensitivity of the model's 
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solution to topography is low in this region. All sensitivity modes correspond- 
ing to small eigenvalues are concentrated in this region. That means, the data 
assimilation procedure intended to reconstruct the topography in this region 
may not be efficient because the influence of corresponding singular modes on 
the model's solution is small. One should use additional apriori information to 
reconstruct topography. 

Turning attention to more realistic models, a number of problems can be 
encountered. First of all this concerns a multi-layer model with different geom- 
etry in each layer. The presence of baroclinic component may also change the 
sensitivity of the model. Second, this study gives no information about par- 
ticular schemes of parametrisation of topography. Numerous modern schemes 
like partial step or shaved cells can not be distinguished in this paper as well as 
different grids like Arakawa's ones. Working with the barotropic vorticity equa- 
tion we can not pay an attention to, for example, C-grid, which is frequently 
used in reahstic models. 

Acknowledgments. All the contour pictures have been prepared by the Grid 
Analysis and Display System (GrADS) developed in the Center for Ocean-Land- 
Atmosphere Interactions, Department of Meteorology, University of Maryland. 

5 Appendix: Numerical resolution 

In order to look for a weak solution of the problem (jlip we perform its varia- 
tional formulation: 

<-^,ip{x,y)>=<A{i;,u;)SLJ,ip{x,y)> + <B{i:,uj) — ,ip{x,y)> (26) 

for any function (p{x, y) € i/p (SI). Here, Hq{U) denotes the linear space of func- 
tions that the square is integrable as well as the square of their first derivatives. 
Functions in this space must vanish on the boundary of the domain. Brackets 
<., .> denote the L2 scalar product: 



The variational formulation ((26|l of the problem (fTTjl allows us to search 
the solution by the finite element method (FEM). 

So far the model (O under consideration is similar to the barotropic one and 
the solution produced by the barotropic model of the North Atlantic typically 
includes a western boundary layer with intense velocity gradients, the advantage 
of refining the triangulation along the western boundary of the domain is rather 
clear. This helps to keep the quality of explicit eddy resolution by the model 
while working with lower number of grid nodes. The comparison of finite ele- 
ments (FE) and finite difference (FD) models performed in [3] revealed that the 
difference arose between simulations by FE and FD techniques can be judged 
as insignificant when the number of FE nodes is about 6 times lower than the 
number of FD ones. 

In spite of the fact that the number of operations per time step and grid 
node is much higher for FE model, the possibility of reducing the number of 
grid points considerably diminish the computational cost of a model run. The 
possibility to have a better precision working with a lower number of grid points 




(27) 
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is very valuable in this work allowing us to perform more detailed study of the 
sensitivity. 

The package MODULEF [3] has been used to perform a triangulation of 
a domain. This package produces quasi-regular triangulation of the domain 
basing on the prescribed grid nodes on its boundary. We require the refining of 
the triangulation near the western boundary and especially in the middle of the 
domain where velocity gradients are extremely sharps. 

The domain is covered by a set of non-intersecting triangles. The set of 
integration points is defined as the union of vertices and mi-edges of triangles. 
Finite elements of type P2 are used here, i.e. the polynomials of the second 
degree Pi{x, y) — aix'^ + bixy + Ciy^ + diX + eiy + fi. The ith finite element is 
taken to be equal to 1 at the i-th integration point and zero at all other points. 

According to the Dirichlet boundary conditions (|14p . we consider internal 
points of the domain only: {xi,yi) G ^l\dfl for i = 1, . . . , N., so variables of the 
problem are presented as linear combinations 

N 



ip{x,y,t) = ^ifj^{t)p^{x,y), 



uj{x, y,t) + /q + Py _^(uj + fQ + I3y 



H{x,y) 



V) i=i V / i j=i 

To simplify notations, we define matrices of mass and rigidity as 

( i = 1 N 

Mi^j = < p^,pj >, C^J = < S7p„\/pj > ^ _ 1 ' " ' ' w (29) 



J =1,...,N 

We multiply the equation pTjl by finite elements pk{x,y), \/k — 1 . . . N 

N 



<qI^ ^^iWPi (a;, y) , Pk {x; y) >=< uj)5uj,pk {x, y)> + 

i—l 

6H 

<B{^,ij) — ,Pk{x,y)> (30) 

Scalar products with operators A and B are developed as follows. From (fT^ 
we get 

<A(v,^)5^,Pfe> = - < J(V, ^),Pfe> + <J{ '^^ i^li^) ^^^'^"^ 

— 1/ <V(5aj, Vpfe> — (T <(5tJ,pfc> 
using (|28|) we write the first Jacobian of (f3T|) as 

bio ^ ^ ( buj\ 

<^(^,^),Pfe>-EE^'(*)(:^) <J{v^.Pi).Vk> (32) 
i=i j=i ^ ' i 
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and the second one as 

N N 



\ J i^l j^l \ J i 

(33) 



where ^ is determined as the solution of the equation 

H 

To solve this equation, we multiply it by Vfc = 1 . . . iV and integrate 

by parts. That gives 

- <— V^, Vpfc>=<(5w,pfe> 
u 

Supposing ^ to be discretised in the same way as other functions = 

AT 

ii[^)Vt{x,y), we get 

N N ^ ^ N 

51 XI ( ) ^1 <P^^P3 ' "^Pk >=^<Pi. Pj > SUJJ 
i=l j = l ^ i j"=l 

or in matricial form 

^ / 1 \ 

m = M5L0 where ^ = XI ( ^7 ) <Pk"^P3^ > (34) 

k=i ^ ^ k 

Combining ^ and ^ we get 

<A{i},uj)5Lu,pk> = {A^^^ ^ A'^^'^'H-^M-vC-aM)5Lu (35) 



i=l 

AT 



- E(^^^^4^) <^(p.P.),P.> (36) 

The scalar products with the operator B in ([50]) is developed in a similar 
way. From (fO|) we get 

r., , , UJ + fo + PySH ^ 

<B{i^,uj)—,pk> = <j(v^, — jr^ir^'P^^~^ 

, <.(^^,(vlv)-'(vlfv.),^., 

Using (P5| we get for the first Jacobian 

<Jw, 77 -^)'Pfe>= z^z^V'»(0( 77 I (-^1 <J{pi.P3)^Pk> 

(38) 
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and for the second one 

<Ji '^)'P'=>= Z^Z^( ) Vj <JiPi,Pj),Pk> (39) 

i=i j=i ^ ^ » 

where 77 is determined from the equation 

vlv), . (vlf V. 

multiplying this equation hy pk{x,y), \fk — 1 . . . N and integrating by parts we 
get 

< j^^V, ^Pk >=< ffjf^'^^ ^Pk > 
The function rj is developed as liner combination of finite elements ri[x,y,t) — 

N 

J2 'ni{t)pi{x,y) and we get 

i=l 

'^j- <p^^pj^^p'^>=J2J2(j^) (^) ^1 <p^^PJ,^Pk> 

i=l j = l ^ ^ i i=l j = l ^ ^ * ^ ^ * 

or in matricial form 

Tirj = V— where Vk,i =^ijj] ^3 <pNPj,^Pk > 
Combining dS?]), (jSHl), dSi) and (001) we get 



(40) 



<S(^,^)^,Pfc> = 5(2)7^-17')^ (41) 

- ^:Mt)(^^^^) <j{p.,P3),Pk> 



i=l 
N 



(2) / uj + fo + f3y\ 

^Ij = 77 ) <J{P^,PJ),Pk> 

i=l ^ / i 

Resuming all of the above, we get the finite element approximation of the 
equation (fTTj) 

^dSi^t) ^ (^A^i) +A^'^H-'M-i^C-aM)6Lo + {B^'^ +b('^H-WS) 
at H 

N 



^k] = E^'W(^) <J{P^,P3),Pk> 



AT 



/i(2) ^ + /o + \ ^ ^/ N ^ 

= j <J{P^,P3),Pk> 
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